2. Métodos Numéricos (2/2)

Autor: Claudio Canales Donoso
Página: ccdonoso.github.io
Cursos: Mecánica Computacional - Diseño Computarizado
Universidad de Santiago de Chile
Fecha: 05/03/24
License: BSD 3 clause
Contenido : Métodos Númericos - 2. - Gradientes Descendientes - Gradientes Conjugados - Sistemas de Ecuaciones Lineales - Cuadratura Gaussiana - Newton Raphson Generalizado
2.1. Gradientes Descendientes
El Descenso de Gradientes [1] es un algoritmo de optimización genérico capaz de encontrar soluciones óptimas a una amplia gama de problemas. La idea general del Descenso de Gradientes es ajustar los parámetros de forma iterativa para minimizar una función de coste \(J\).
Supongamos que usted se encuentra perdido en las montañas en medio de una densa niebla, y que sólo puede sentir la pendiente del suelo bajo sus pies. Una buena estrategia para llegar rápidamente al fondo del valle es ir cuesta abajo en la dirección de la pendiente más pronunciada. Esto es exactamente lo que hace el Descenso del Gradiente: mide el gradiente local de la función de error con respecto al vector de parámetros \(θ\), y va en la dirección del gradiente descendente. Una vez que el gradiente es cero, ¡se ha alcanzado un mínimo!
Concretamente, se empieza llenando \(\boldsymbol{\theta}\) con valores aleatorios (esto se llama inicialización aleatoria). A continuación, se mejora gradualmente, dando un pequeño paso a la vez, cada paso tratando de disminuir la función de coste (por ejemplo, el MSE), hasta que el algoritmo converge hasta que el algoritmo converja a un mínimo .

El gradiente siempre esta orientado hacia la dirección de máximo crecimiento
El gradiente de una función siempre esta orientado a la dirección de máximo crecimiento y la dirección opuesta es la dirección de mínimo crecimiento. En base a esta propiedad del gradiente, se puede utilizar la dirección y magnitud del gradiente para encontrar un máximo o un mínimo.
El método consiste en los siguientes pasos:
Inicialización de un vector de parámetros \(\boldsymbol{\theta}_0\)
Repetir hasta que un criterio de término se cumpla:
\(\boldsymbol{\theta}_{n+1} = \boldsymbol{\theta}_n - \alpha \triangledown J(\boldsymbol{\theta}_n), \quad n>0\)
Donde \(\alpha\) pertenece a \(\mathbb{R}^{+}\) y tiene que ser lo suficientemente pequeño para que cumpla que \(J(\boldsymbol{\theta}_{n+1}) < J(\boldsymbol{\theta}_n)\)
En Resumen el método se basa en descender de la cumbre utilizando los gradientes de la función dando pequeños pasos.
Recordar que este método esta condicionado a la inicialización.
2.2. Resolver el siguiente problema de optimización
[5]:
%matplotlib inline
from IPython.display import HTML
import matplotlib.pyplot as plt
import matplotlib.animation
import numpy as np
# Datos de la Grafica de la función
x = np.linspace(0,8) #ndarray
f = lambda x: x**2. * np.sin(x)
y = f(x)
# Setup Graficos
fig, ax = plt.subplots(figsize=(8,6))
ax.axis([0,8,-30,30])
l, = ax.plot([],[],'rx')
time_text = ax.text(3, 20, '', fontsize=14)
ax.plot(x,y)
# Ingresa el Código aquí
# ---------------------------
#Gradiente Descendiente
learning_rate = .001
grad = lambda x: 2.*x*np.sin(x)+x**2.*np.cos(x)
theta_n = [np.pi]
iters = 200
for i in range(iters):
theta_n.append(theta_n[i] - learning_rate * grad(theta_n[i]))
theta_n = np.array(theta_n)
J = f(theta_n)
#-------------------------------
def animate(i):
time_text.set_text(r"El valor de $\theta$: {:3f}".format(theta_n[i])+" \nIteración: {}\nFunción objetivo:{:2f}".format(i,J[i]))
l.set_data(theta_n[:i], J[:i])
return time_text
ani = matplotlib.animation.FuncAnimation(fig, animate, frames=len(J))
ax.legend(["Iteraciones del método","Función a minimizar"])
plt.close()
HTML(ani.to_jshtml())
[5]:
2.2.1. Resolver el siguiente problema con gradientes descendientes.
Hacer un código, que sea capaz de resolver el siguiente problema de minimización para cualquier \(n\) y con inicialización \(\mathbf{0}\).
Para el caso \(n=2\) se tiene que la función es:

[6]:
def sum_x2(x,n):
"""
Esta función calcula la función a minimizar
Parameter
----------
x : ndarray
El argumento de la función
n : int
Determina la dimensión del problema
Return
----------
f: float
El valor de la función
"""
f = 0.
for i in range(n):
f = (x[i]-1.)**2 + f
return f*0.5
def grad_sum_x2(x,n):
"""
Esta función calcular el gradiente
Parameter
----------
x : ndarray
El argumento de la función
n : int
Determina la dimensión del problema
Return
----------
grad: ndarray
"""
grad = np.zeros(n)
for i in range(n):
grad[i] = (x[i]-1.)
return grad
n = 40000
iterations = 100
learning_rate = 0.1
x = np.zeros(n)
for i in range(iterations):
x = x - learning_rate * grad_sum_x2(x,n)
print(x)
[0.99997344 0.99997344 0.99997344 ... 0.99997344 0.99997344 0.99997344]
2.2.2. Solución de Sistemas de ecuaciones Lineales a través de optimización
Un sistema de ecuaciones puede ser escrito como
,
donde \(\mathbf{A}\) es una matriz \(n\times n\), las expresiones \(\mathbf{x}\) y \(\mathbf{b}\) son matrices de \(n\times 1\).
Este problema puede ser escrito de la siguiente manera:
La existencia de un minimizador esta garantizada si es que \(\mathbf{A}\) es una matriz simetrica definida positiva.
2.3. Método de los Gradientes Conjugados
En matemáticas, el método del gradiente conjugado es un algoritmo para la solución numérica de sistemas de ecuaciones lineales determinados, es decir, aquellos cuya matriz sea positiva-definida. El método del gradiente conjugado suele implementarse como un algoritmo iterativo, aplicable a sistemas dispersos que son demasiado grandes para ser tratados por una implementación directa u otros métodos directos como la descomposición de Cholesky. Los sistemas dispersos de gran tamaño suelen surgir al resolver numéricamente ecuaciones diferenciales parciales o problemas de optimización.
El método del gradiente conjugado también puede utilizarse para resolver problemas de optimización sin restricciones, como la minimización de la energía.
El método del gradiente biconjugado proporciona una generalización a las matrices no simétricas. Varios métodos de gradiente conjugado no lineal buscan los mínimos de las ecuaciones no lineales y las funciones objetivo de caja negra.
Para poder utilizar este metodo, la matriz tiene que ser definida positiva
Para asegurar que una función presenta un mínimo local en \(\mathbf{x}\), se tiene que satisfacer que
Es fácil comprobar que todas las condiciones se cumplen para un sistema de ecuaciones con una matriz simétrica y definida positiva.
Se observa que al encontrar el valor que minimiza la función, al mismo tiempo se resuelve el sistema de ecuaciones lineales y la condición de optimalidad se verifica con la condición del Hessiano.
2.3.1. Este método utiliza los gradientes descendientes para encontrar el óptimo del problema, pero los pasos que da son conjugados de su paso anterior ( Gradientes Conjugados).
Como se puede ver, el gradiente esta definido y el primer vector de la base es \(\mathbf{p_0}\) que toma el valor negativo del gradiente en el primer paso \(\mathbf{p_0} = \mathbf{b-Ax}\). Los otros vectores de la base serán conjugados de este. Notar que este valor tambien es el residuo del método.
Como se puede observar de las ecuaciones anteriores, \(\mathbf {r} _{n}\) es el valor negativo del gradiente de \(f\) en $ \mathbf {x} {k}$, por lo que el método de gradiente descendente requeriría moverse en la dirección :math:`mathbf{r}_k`. Aquí, sin embargo, insistimos en que las direcciones $ :nbsphinx-math:`mathbf {p}`{n}$ sean conjugados entre sí. Una forma práctica de aplicar esto es exigir que la siguiente dirección de búsqueda se construya a partir del residuo actual y de todas las direcciones de búsqueda anteriores. La restricción de conjugación es una restricción de tipo ortonormal y, por tanto, el algoritmo puede considerarse un ejemplo de ortonormalización de Gram-Schmidt. Esto da la siguiente expresión:
Siguiendo esta dirección, la siguiente ubicación óptima viene dada por:
con
2.3.2. El algoritmo para resolver un sistema de ecuaciones lineales siendo \(\mathbf{A}\) definida positiva es:
Es posible observar que este algoritmo es una versión modificada del método de los gradientes descendientes
2.4. Cuadratura Gaussiana
En el análisis numérico, una regla de cuadratura es una aproximación de la integral definida de una función, normalmente expresada como una suma ponderada de los valores de la función en puntos específicos dentro del dominio de integración. Una regla de cuadratura gaussiana de n puntos, llamada así en honor a Carl Friedrich Gauss, es una regla de cuadratura construida para obtener un resultado exacto para polinomios de grado \(2n - 1\) o menos mediante una elección adecuada de los nodos \(x_i\) y los pesos \(w_i\) para \(i = 1,\dots ,n\). La formulación moderna que utiliza polinomios ortogonales fue desarrollada por Carl Gustav Jacobi en 1826 .El dominio de integración más común para dicha regla se toma como \([-1, 1]\), por lo que la regla se establece como:
Para el problema de integración más simple indicado anteriormente, es decir, \(f(x)\) es bien aproximado por polinomios en \([-1,1]\), los polinomios ortogonales asociados son polinomios de Legendre, denotados por \(P_n(x)\). Con el n-ésimo polinomio normalizado para dar \(P_n(1) = 1\), el i-ésimo nodo de Gauss, \(x_i\), es la i-ésima raíz de \(P_n\). Estos pesos se pueden calcular como:

A través de una transformación lineal, es posible cambiar los límites de integración de la siguiente manera:
Calcular la siguiente integral utilizando Cuadratura Gaussiana :
[7]:
def gauss_quadrature3(f,a,b):
"""
Esta función integra utilizando 3 puntos de cuadratura Gaussiana
Parámetros
----------
f : function
función a integrar
a : float
Límite inferior del intervalo de integración
b : float
Límite superior del intervalo de integración
Retorna
----------
integral : float
Valor de la approximación de la integral
"""
wi = np.array([5./9.,8./9.,5./9.])
xi = np.array([-np.sqrt(3./5.),0.,np.sqrt(3./5.)])
integral = 0.
for i in range(3):
integral = wi[i]*f(xi[i]*(b-a)*0.5+(a+b)*0.5) + integral
return integral*(b-a)*0.5
print(gauss_quadrature3(lambda x: x**3. + x**5. + np.sin(x),-2.,1.))
-15.206934548042273
2.5. Desafío: Programar la cuadratura Gaussiana para cualquier \(n\).
[8]:
def gauss_quadrature(f,a,b,n=3):
"""
Esta función integra utilizando n puntos de cuadratura Gaussiana
Parámetros
----------
f : function
función a integrar
a : float
Límite inferior del intervalo de integración
b : float
Límite superior del intervalo de integración
n : int
Cantidad de puntos de la cuadratura Gaussiana
Retorna
----------
integral : float
Valor de la approximación de la integral
"""
2.6. Método de Newton Raphson Generalizado
El método de Newton Raphson puede ser utilizado para solucionar sistemas de ecuaciones no-lineales.
donde
Tambien se puede definir un vector de funciones
Por lo tanto el sistema de ecuaciones puede ser escrito como:
Las series de Taylor nos permiten expresar cada una de las funciones como:
Es posible utilizar el termino de primer orden para linealizar la ecuación, es decir hasta el gradiente.

Donde:
Si se dan cuenta, ahora el sistema de ecuación es un lineal y \(\mathbf{f(x)=0}\). Por lo tanto, al linealizar el problema, es posible realizar pequeñas aproximaciones en base a vector de inicialización \(\mathbf{x}_0\) y iterar sucesivamente al igual que el método de Newton-Raphson 1D. Por lo tanto, con un vector \(\mathbf{x}_0\) se puede obtener un vector \(\mathbf{x}_1=\mathbf{x}_0 + \mathbf{\Delta x_0}\) y ese vector se puede obtener al solucionar el siguiente sistemas de ecuaciones:
El algoritmo básico consiste en repetir estos tres pasos:
Resolver el sistema de ecuaciones:
Esta ecuación es semejante a resolver:
Actualizar el vector solución:
Repetir estos pasos hasta que un criterio de término se cumpla:
Revisar este documento que contiene ejemplos de como utilizar el método
Resolver el siguiente problema con \(\mathbf{x_0}=(0,0)\)
[9]:
# Hacer el código aquí
def vector_funciones(x):
vf = np.zeros(2)
vf[0] = 5*x[0]**2+3*x[0]*x[1]-2
vf[1] = x[0]**2+7*x[1]**2+3*x[0]*x[1]-10
return vf
def jacob(x):
"""
df1/dx df1/dy
df2/dx df2/dy
"""
J = np.zeros((2,2))
J[0,0] = 10*x[0] + 3*x[1]
J[0,1] = 3*x[0]
J[1,0] = 2*x[0]+3*x[1]
J[1,1] = 14.*x[1]+3.*x[0]
return J
x = np.array([-0.1,-0.1])
for i in range(20):
dx = np.linalg.solve(jacob(x), -1.*vector_funciones(x))
x = dx + x #xn+1
residuo = vector_funciones(x)
norma = np.dot(residuo,residuo)
print("El valor del residuo:",norma)
El valor del residuo: 3.1554436208840472e-30
2.7. Problema : encontrar la posición de un avión utilizando tres antenas y el método de Newton Raphson
Para monitorear la posición de un avión durante el vuelo, este se comunica de forma simultánea con tres antenas, de forma tal que de cada una de ellas obtiene el tiempo que tarda una señal en ir desde el avión a las antenas. Considerando que las señales viajan a \(299.792 km/s\), y utilizando el archivo “tiempo_avion.dat”, encuentre la posición instantánea \(x,y,z\) del avión, para cada punto de muestreo, utilizando el método de newton Raphson en varias variables.

Se le pide guardar la posición \(x,y,z\) del avión en un archivo de texto llamado “posicion.dat”
Las distancias se pueden calcular como: \(d=ct\), donde \(c\) es la velocidad de la luz y \(t\) es el tiempo de ida de la señal.
La posición de las antenas en [km] son:
Consejo, utilice como punto de inicialización :math:`(0,0,0)` :math:`[km]`
[10]:
import numpy as np
def vector_funciones(x,d):
vf = np.zeros(3)
x , y , z = x
d1 , d2, d3 = d
vf[0] = x**2+(y-50.)**2+z**2.-d1**2.
vf[1] = (x-60.)**2 +y**2 + z**2 -d2**2
vf[2] = (x-70.)**2 + (y-30.)**2 + z**2-d3**2.
return vf
def jacob(x,d):
"""
df1/dx df1/dy df1/dz
df2/dx df2/dy df2/dz
df3/dx df3/dy df3/dz
"""
x , y , z = x
d1 , d2, d3 = d
J = np.zeros((3,3))
J[0,0] = 2*x
J[0,1] = 2.*(y-50.)
J[0,2] = 2.*z
J[1,0] = 2.*(x-60.)
J[1,1] = 2.*y
J[1,2] = 2.*z
J[2,0] = 2.*(x-70.)
J[2,1] = 2.*(y-30.)
J[2,2] = 2*z
return J
tiempos_s = np.loadtxt(fname="files/tiempo_avion.dat")/1000. # columnas están ordenadas como[t1,t2,t3]
# Esta expresión se multiplica por 1000
# para que el tiempo este en Segundos.
vel_luz = 299792. # Vel Luz km/s
distancias = tiempos_s*vel_luz # "Tiempo_señal x Velocidad"
mediciones = tiempos_s.shape[0] # int -> Cantidad de filas -> Cantidad de mediciones de los radares
x = np.array([0.1,0.1,0.1])
sol = np.zeros((mediciones,3))
for t in range(mediciones):
for i in range(20):
dx = np.linalg.solve(jacob(x,distancias[t,:]), -1.*vector_funciones(x,distancias[t,:]))
x = dx + x #xn+1
sol[t,:] = x # observar que el x obtenido del paso anterior
# se utiliza como semilla para t+1
2.8. Solución numérica con la información del avión.
[11]:
import matplotlib.pyplot as plt
import mpl_toolkits
from mpl_toolkits.mplot3d import Axes3D
from numpy.random import rand
from IPython.display import HTML
from matplotlib import animation
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111, projection='3d')
xline , yline, zline = sol.T
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
def animate(frame):
ax.view_init(30, frame)
ax.plot3D(xline[:frame], yline[:frame], zline[:frame], color = 'green')
return fig
ax.legend(["Trayectoria de avión"])
plt.close()
anim = animation.FuncAnimation(fig, animate, frames=len(zline))
HTML(anim.to_jshtml())
[11]:
2.8.1. Solución analítica
[12]:
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from numpy.random import rand
from IPython.display import HTML
from matplotlib import animation
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111, projection='3d')
zline = np.linspace(0, 15, 100)
xline = np.sin(zline)
yline = np.cos(zline)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
def animate(frame):
ax.view_init(30, frame)
ax.plot3D(xline[:frame], yline[:frame], zline[:frame], color = 'blue')
return fig
ax.legend(["Trayectoria de avión"])
plt.close()
anim = animation.FuncAnimation(fig, animate, frames=len(zline))
HTML(anim.to_jshtml())
[12]:
2.9. Referencias
[1] Géron, A. (2019). Hands-on machine learning with Scikit-Learn, Keras, and TensorFlow: Concepts, tools, and techniques to build intelligent systems. O’Reilly Media.